#setwd("filepath")
path<-getwd()
library(foreign)
library(rjags)
set.seed(666)

## Read in IACRT acceptance data
iacrt.dat<-read.table(file="iacrt.csv",sep=",",header=T)

## Read in ICCPR ratification data
ccpr.dat<-read.table(file="iccpr.csv",sep=",",header=T)
ccpr.dat<-ccpr.dat[ccpr.dat$cowcode<200&ccpr.dat$year<2005,]

## Read in latent democracy data
dem.dat<-read.csv("uds_summary.csv", header=TRUE, sep=",")
#lji.dat<-read.csv("LJI-estimates.csv", header=TRUE, sep=",")

## Merge

new.dat<-merge(iacrt.dat, ccpr.dat, by.x=c("ccode", "year"), by.y=c("cowcode", "year"), all.x=T)

new.dat<-merge(new.dat, dem.dat, by.x=c("ccode", "year"), by.y=c("cowcode", "year"), all.x=T, all.y=T)
new.dat<-new.dat[new.dat$ccode<200&new.dat$year>1960,]

new.dat<-new.dat[order(new.dat$ccode, new.dat$year),]

## Mark first observation in each panel
new.dat$new.panel<-NA
for (i in 6:dim(new.dat)[1]){
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-1], 1, 0)
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-2], 1, 0)
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-3], 1, 0)
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-4], 1, 0)
	new.dat$new.panel[i]<-ifelse(new.dat$ccode[i]!=new.dat$ccode[i-5], 1, 0)
}
new.dat$new.panel[1:5]<-0

post.draws<-matrix(nrow=dim(new.dat)[1],ncol=10)
for (i in 1:dim(new.dat)[1]){
	for (j in 1:10){
		post.draws[i,j]<-rnorm(1, mean=new.dat$mean[i], sd=new.dat$sd[i])		
	}
}

## Difference in latent democracy
difs<-matrix(nrow=dim(new.dat)[1],ncol=50)
for (i in 6:dim(post.draws)[1]){
	for (j in 1:10){
		difs[i,j]<-post.draws[i,j]-post.draws[i-1,j]
		difs[i,j+10]<-post.draws[i,j]-post.draws[i-2,j]
		difs[i,j+20]<-post.draws[i,j]-post.draws[i-3,j]
		difs[i,j+30]<-post.draws[i,j]-post.draws[i-4,j]
		difs[i,j+40]<-post.draws[i,j]-post.draws[i-5,j]
	}
}

new.dat<-cbind(new.dat,post.draws,difs)

names(new.dat)[16:75]<-c("draw.1","draw.2","draw.3","draw.4","draw.5",
"draw.6","draw.7","draw.8","draw.9","draw.10",
"draw.1.dif","draw.2.dif","draw.3.dif","draw.4.dif","draw.5.dif",
"draw.6.dif","draw.7.dif","draw.8.dif","draw.9.dif","draw.10.dif",
"draw.1.dif.2","draw.2.dif.2","draw.3.dif.2","draw.4.dif.2","draw.5.dif.2",
"draw.6.dif.2","draw.7.dif.2","draw.8.dif.2","draw.9.dif.2","draw.10.dif.2",
"draw.1.dif.3","draw.2.dif.3","draw.3.dif.3","draw.4.dif.3","draw.5.dif.3",
"draw.6.dif.3","draw.7.dif.3","draw.8.dif.3","draw.9.dif.3","draw.10.dif.3",
"draw.1.dif.4","draw.2.dif.4","draw.3.dif.4","draw.4.dif.4","draw.5.dif.4",
"draw.6.dif.4","draw.7.dif.4","draw.8.dif.4","draw.9.dif.4","draw.10.dif.4",
"draw.1.dif.5","draw.2.dif.5","draw.3.dif.5","draw.4.dif.5","draw.5.dif.5",
"draw.6.dif.5","draw.7.dif.5","draw.8.dif.5","draw.9.dif.5","draw.10.dif.5")

for (i in 26:75){
	new.dat[,i][new.dat$new.panel==1]<-NA
	}

new.dat<-new.dat[new.dat$year>1965,]

civwar.data<-read.csv(file="ucdp_onset2012.csv", sep=",", header=T)
cw.dat<-data.frame(cbind(civwar.data$year, civwar.data$gwno, civwar.data$incidencev412))
colnames(cw.dat)<-c("year", "ccode", "conflict")
cw.dat<-cw.dat[order(cw.dat$ccode, cw.dat$year),]
cw.dat<-cw.dat[cw.dat$ccode<200&cw.dat$ccode>0&cw.dat$year<=2004,]

newer.dat<-merge(new.dat, cw.dat, by=c("ccode", "year"), all.x=T)

legal.data<-read.dta(file="mitchell_jpr.dta")
legal.dat<-legal.data[legal.data$cowcode<200&legal.data$year>=1965&legal.data$year<=2004,]
legal.dat<-with(legal.dat, cbind(cowcode, year, common))

newest.dat<-merge(newer.dat, legal.dat, by.x=c("ccode", "year"), by.y=c("cowcode","year"),all.x=T)

## fill in missing legal system data
newest.dat$common<-ifelse(is.na(newest.dat$common), 0, newest.dat$common)
common<-unique(newest.dat$ccode[newest.dat$common==1])
for (i in 1:length(common)){
	newest.dat$common<-ifelse(newest.dat$ccode==common[i], 1, newest.dat$common)
}

newest.dat<-newest.dat[newest.dat$year<2005,]

save(newest.dat,file="ccpr_data_NAs.rda")

dat<-newest.dat[,-5]
dat<-dat[,-5]
dat<-na.omit(dat)
dim(dat)

## Create time counter

dat$end<-ifelse(dat$year==2004,1,0)
dat$end<-ifelse(dat$iccpr_rat==1,1,dat$end)
dat$id<-rev(cumsum(rev(dat$end)))

dat<-dat[order(dat$id, dat$year),]
helper<-rep(1,dim(dat)[1])
dat$t<-as.vector(unlist(tapply(helper,dat$id,cumsum)))
#dat$t2<-dat$t^2
#dat$t3<-dat$t^3

common.l2<-as.numeric(tapply(dat$common, dat$id, mean))

## 442 country-years 
nrow(dat)

### JAGS data and model

the.hi.model.list<-list()
for (i in 1:10){

	N<-nrow(dat)
	id<-dat$id
	ccpr<-dat$iccpr_rat
	dem<-dat[,i+13]
	conflict<-dat$conflict
	common<-common.l2
	t<-dat$t
#	t2<-dat$t2
#	t3<-dat$t3

	the.data<-list("N"=N, "id"=id, "ccpr"=ccpr, "dem"=dem, "conflict"=conflict, "common"=common, "t"=t)

	the.inits.1<-list(b=c(0, 0, 0, 0), a=rep(0,28), tau=1, .RNG.seed=666, .RNG.name="base::Super-Duper")
	the.inits.2<-list(b=c(1, -1, -1, 1),a=rep(-1,28), tau=5, .RNG.seed=666, .RNG.name="base::Super-Duper")
	the.inits<-list(the.inits.1,the.inits.2)

	the.hi.model<-jags.model("ccpr_hi_logit.bug", data=the.data, n.chains=2, n.adapt=0, 	inits=the.inits)
	update(the.hi.model,10000)
	the.hi.model.list[[i]]<-coda.samples(the.hi.model, c("a", "b", "g"), 10000, thin=1)
}

save(the.hi.model.list, file="ccpr_demsq_models.rda")

the.hi.model.dif.list<-list()
for (i in 1:50){

	N<-nrow(dat)
	id<-dat$id
	ccpr<-dat$iccpr_rat
	dif<-dat[,i+23]
	conflict<-dat$conflict
	common<-common.l2
	t<-dat$t
#	t2<-dat$t2
#	t3<-dat$t3

	the.data<-list("N"=N, "id"=id, "ccpr"=ccpr, "dif"=dif, "conflict"=conflict, "common"=common, "t"=t)

	the.inits.1<-list(b=c(0, 0, 0), a=rep(0,28), tau=1, .RNG.seed=666, .RNG.name="base::Super-Duper")
	the.inits.2<-list(b=c(1, -1, 0.01),a=rep(-1,28), tau=5, .RNG.seed=666, .RNG.name="base::Super-Duper")
	the.inits<-list(the.inits.1,the.inits.2)

	the.hi.model<-jags.model("ccpr_hi_logit_dif.bug", data=the.data, n.chains=2, n.adapt=0, inits=the.inits)
	update(the.hi.model,20000)
	the.hi.model.dif.list[[i]]<-coda.samples(the.hi.model, c("a", "b"), 10000, thin=1)
}

save(the.hi.model.dif.list, file="ccpr_dif_models.rda")


## Summaries, check convergence
lapply(the.hi.model.list,summary)
lapply(the.hi.model.list,gelman.diag)

lapply(the.hi.model.dif.list,summary)
lapply(the.hi.model.dif.list,gelman.diag)

